en 



> 

X 



Measuring Directionality in Random Fields 
Using the Monogenic Signal 

Sofia Olhede, David Ramirez, Member, IEEE and Peter J. Schreier, Senior Member, IEEE 



Abstract 



^^ Analyzing directional structures in images is important in many applications since one-dimensional patterns often 

correspond to important features such as object contours or trajectories. Most of the work on this topic has dealt 

O , with deterministic images; so far, little research has been performed on analyzing directionality in random fields. In 

'^^ this paper, we propose a measure to quantify the degree of directionality in random fields. This measure can be used 

("^ to construct statistical tests for the presence of directional structure. Our measure is based on the monogenic signal, 

which enables a unique decomposition of a 2D signal into local amplitude, local orientation, and local phase. We 

I I investigate the second-order statistical properties of the monogenic signal for isotropic, anisotropic, and unidirectional 

V \ random fields. We analyze our measure of directionality for finite-size sample images, and determine a threshold to 

t/3 distinguish between unidirectional and non-unidirectional random fields. 

" Index Terms 

> 

00 

^v Anisotropy, monogenic signal, quaternions, Riesz transform, stationary random field, unidirectional. 

.• I. Introduction 

^^ This paper introduces a novel measure of the degree of directionality of a random field. The analysis of directional 

structure in images is crucial to many applications since one-dimensional patterns often correspond to important 
features such as object contours or trajectories. Identifying the presence of anisotropy and estimating the orientation 
of such patterns is particularly important; see, for instance, the detection of ship wakes in Synthetic Aperture Radar 
(SAR) images |1|, the analysis of texture and optical flow by estimating multidimensional orientation p), the 
efficient coding of local differential structures in images |]3l, and the analysis of superimposed directional patterns, 
which may occur in X-ray projection imaging |4). However, all the techniques developed in the aforementioned 
papers are deterministic, and there exists very little work on analyzing directional structures in random fields. 

A random field is a stochastic process whose argument is a multidimensional vector In our case, the argument is a 
2D vector, and a spatial 2D random field may also be called a random image. Random fields have ubiquitous practical 
applications in areas as diverse as geophysics 15), Q, oceanography 17), and medical imaging 181. A fundamental 
characteristic of a spatial random field is the degree of rotational invariance of its second-order statistical properties 

This work was presented in part at the 2012 IEEE International Conference on Image Processing (ICIP). 
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||9] p. 57]. A spatially isotropic random field exhibits perfect rotational invariance as its spatial covariance displays 
circular contour levels. Fields that are not isotropic are called anisotropic. A common subclass is the class of 
geometrically anisotropic fields |[9] p. 61] that have covariances with elliptical rather than circular contours. Taking 
this idea to the extreme, we arrive at unidirectional random fields, where there exists a rotation so that there is 
variation only in one of the two axes. Figure [T] shows samples of three random images: isotropic, geometrically 
anisotropic, and unidirectional, giving insight into the features present in a field when departing from the assumption 
of isotropy. 

A large fraction of work in statistics has focussed on isotropic random fields — which is at least partially due 
to the fact that these have convenient mathematical properties — and many spatial models utilize isotropy ilO|. Yet 
it is also clear that there are many examples where isotropy is a very poor assumption: For instance, wind data 
or ocean currents can be highly directional. In order to deal with anisotropic random fields, there are a number 
of models ||9l-| 1 1 1, the most common of which is geometric anisotropy. However, for highly directional random 



fields, it would seem that the most appropriate model should be one that consists of unidirectional components. 

In the deterministic case, there exist quite a few papers that deal with the estimation of local orientation, e.g., Q- 
p), p2) . For a brief summary of some elementary techniques, we refer the reader to the background section in | [T2) . 
We use a different approach based on the random monogenic signal, which allows us to define a statistical measure 
for the degree of unidirectionality and to construct statistical tests for the presence of directional structure. The 



monogenic signal 1 13 enables the unique decomposition of a two-dimensional real image /(x), with x = [xi, X2] , 



into a local amplitude, a local orientation, and a local phase. It is arguably the most compelling 2D-generalization 
of the analytic signal, which enables the unique decomposition of a real signal x{t) — a{t) cos (f>{t) into local 
amplitude a{t) and local phase 0(i). The monogenic signal has received considerable attention, with applications in 
image processing ranging from contour detection and local structure analysis p4) , fTS) , stereo, motion estimation, 
and image registration p6) to image segmentation and phase-contrast imaging fVJ]. There have also been extensions 



of the concept to a multiresolution monogenic signal in the wavelet domain p8) , 1 19|. 

As we will review in Section |ll] the monogenic signal is constructed by complementing the original signal 
/(x) with its two Riesz transforms g(x) and /i(x) pO[ . It can either be represented as a three-dimensional vector 
[f{x),g{x), h{x)]'^ or as a quaternion m(x) = /(x) + i.g(x) + jh{x.) + fc • 0, where the fc-part remains zero. So 
far, work on monogenic signals has focussed only on the deterministic case. Aside from our own conference paper 
pT), we are not aware of any work that has been performed on a random monogenic signal. In Section III we 



investigate the second-order statistical properties of random monogenic signals for stationary random fields, and in 



Section IV we examine these properties for the special and important cases of isotropic, geometrically anisotropic, 
and unidirectional random fields. 

In Section W\ we propose and analyze a measure to quantify the degree of unidirectionality for a random field. 
In particular, we show that, for an infinite-size sample image, our measure is identically one only for unidirectional 
random fields and zero for isotropic random fields, with values in between for other degrees of anisotropy. For a 
finite-size sample, our measure of directionality is no longer guaranteed to be one. We determine its finite sample 
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expectation, by carefully expanding the properties depending on the size of the image, and a threshold to distinguish 



between unidirectional and non-unidirectional fields. Finally, in Section VI we illustrate the performance of our 
measure on simulated and real random fields. 

II. Preliminaries 
A. Quaternion random vectors 

In the one-dimensional case, it is common practice to encode the signal and its Hilbert transform in one complex- 
valued analytic signal. While this is not strictly necessary (obviously one could also work with a 2D vector instead), 
this practice is universally accepted because it illuminates and simplifies matter significantly. If we would like to 
do something similar with the monogenic signal, we need to employ the algebra of quaternions. 

In this section, we provide a brief review of quaternion algebra and the second-order analysis of quaternion 
random vectors (see, e.g., p2[). Quaternions are 4D hypercomplex numbers, first proposed by Hamilton |J23), and 
defined as 

q — a -\- bi -\- cj -\- dk, (!) 

where a, b, c, and d are real numbers and i, j, and k are imaginary units satisfying 

ij = k^ -ji jk^i^ -kj (2) 

ki = j = —ik i = j = k = ijk = — 1. (3) 

It is easy to check that quaternions form a normed division algebra H that is non-commutative, i.e., for qi,q2 G H, 
generally qiq2 ^ q2Q_i- The conjugate of q is defined as g* = a — bi — cj — dk and the norm is \q\ = ^/qq* = 
\Jd} + b'^ + c'^ + d^, which satisfies \qiq2\ = |'Zi||'?2|- The inverse of q, with q ^ 0, is q^^ = q* /\q\^, the inner 
product between qi and q2 is Kt{qiq2), and two quaternions are orthogonal if their inner product is zero. Finally, the 
involution of q over a pure unit quaternion -q is given by q^''^ = —VQV- For ^ more complete review of quaternions, 
we refer the reader to ||24|. 

The second-order statistical analysis of a zero-mean complex random vector x is based on the covariance matrix 
R-x.x = cov{x, x} ~ E[x.x-^] and the complementary covariance matrix Rx.x* = cov{x, x*} = i?[xx-'"] 125). 
A common assumption in complex-valued signal processing is propriety, which is characterized by vanishing 
complementary covariance. This can be visualized as rotational invariance because, in the proper case, x and 
xe^" have the same second-order moments for arbitrary real angle a. While propriety can often be justified, there 
are also many situations where it is a very poor model of the underlying physics. 

To completely characterize the second-order statistics of a zero-mean quaternion random vector q, we need the 
covariance matrix and three complementary covariance matrices |l22l. There is some freedom in how to choose 
these complementary covariance matrices. We will employ the most useful choice 

Rq.q = cov{q, q}, Rq,q(,) , = cov{q, q^"*)} Rq_q(V) , = cov{q, q(''')} Rq_q(V') = cov{q, q^''"^}, (4) 
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where 77, rj and 77" are three orthogonal pure unit quaternions, for instance, r\ ^ i, r( = j, rj" = k. Because there 
are three complementary covariance matrices, there are different kinds of propriety for quaternion random vectors 
p2). The only kind of interest to us is C'-propriety. A quaternion random vector is C-proper if and only if both 



R (^/) and R (r,") vanish. This obviously depends on an appropriate choice of rj. If q is C'-proper for some 
?/, it is generally improper for a different choice of 77. In our case, the choice of 77 will be linked to the directional 
structure of the random field. 



B. Monogenic signal 

The monogenic signal was introduced by Felsberg and Sommer in p3) , even though it had already seen some 
prior use in applied mathematics p6| , p7) and geophysics p8) . It is a convenient method of defining an amplitude 
and a vector- valued phase at any point in space |13|, and can be considered as an appropriate generalization to the 
analytic signal pO) . The continuous monogenic signal is defined using the continuous Riesz transform |20|. The 
Riesz transform enjoys a number of convenient mathematical properties, chief among them the commutativity with 
spatial translations and dilations, and equivariance with respect to rotation | pO[ . These properties make the Riesz 
transform a very compelling candidate for separating information about structural and energetic aspects of a 2D 
signal. There are two Riesz transforms in 2D, which araH 

^''V(x) = 7^ ff /(y) /' ~ ""i' dy ^ (rW . /)(x), ; = 1,2, (5) 

27rJJ ||x-y||'' 

where r'^'-'(x) = — a;;/(27r||x||'^), x = [a;i, 2:2]"'", and y = [2/1,2/2]"^- The monogenic signal is formed by placing the 

two Riesz transforms into the i- and j-parts of a quaternion-valued signal, leaving the /c-part empty: 

m(x) = X/(x) = /(x) + ^.g(x) + j/i(x), (6) 

where g{x) = TZ''^' f{x) and /i(x) = TZ^^' f{x). The Riesz transforms are most easily described in the 2D Fourier 
domain, where the Riesz transform kernels are 

R'^^^k) = -i-^ ^ -icos(K) (7) 

i?(2)(k) = -i-^^-isin(K). (8) 

l|k|| 

In this equation, k = [fci,fc2]"^ = fc[cos(K), sin(K)]-'" is the 2D wavenumber, and k = argk is the corresponding 
angle in the 2D plane. If the transformation is implemented over the entire real plane, the monogenic signal can 
be calculated either in the spatial domain using (|5]l or, equivalently, in the 2D Fourier domain using (|7| and ([8]). 

For a sampled lowpass image {fn,n' — f{^n,n'), n, v! — ~N/2, . . . , N/2 — l}rjassumed w.l.o.g. to be sampled 
at a unit sampling period, the Riesz transforms have to be calculated discretely. Analogously to the discrete-time 
Fourier transform (DTFT) and the discrete Fourier transform (DFT), we need to define both a discrete-space Riesz 

'By default, when we do not state the limits of an integral, they are — oo and oo. 
^Without loss of generality, we suppose that A'^ is even. 
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transform, which is continuous in the wavenumber domain, and a discrete Riesz transform, which is discrete and 
periodic in both space and wavenumber domains. To emphasize the difference between the two, we sometimes refer 
to the discrete Riesz transform as the periodic Riesz transform. 

The impulse response of the discrete-space Riesz transform for a non-periodic discrete-space signal is 



r(')(x„,„0 = -7^ r r^#^^^*'"^'^"'^'^k, ? = 1,2. (9) 

{2ttY J_^ J_^ ^kf + kl 

The first Riesz transform of the random field /„.„' is obtained as the convolution 

oo 

<?(x„,„0= Y. ^^'Hxi-„,r-„')/(xi,i')- (10) 

^,Z' — — oo 

and the second Riesz transform /i(x„.„') is obtained analogously by convolving /(x„„/) with r'^^^(x„„/). Using 
the 2D DFT 

N/2-1 

F(k)= Y. /(x„,„')e"'''^"+'^"'\ (11) 

n,n' = -N/2 

with kl E (— TT, tt], the periodic impulse response of the discrete Riesz transform is 

fWfx ,)- — V ^' p>(fci»+fc2»') 7-12 (12) 

i^n.n')- ^2Z.^/^^^^e , ;-l,2. (12) 

The first Riesz transform of /(x„.„/) is then obtained as the circular convolution 

N/2-1 

3(X„,„') = Y '''^^H^l-n mod N, I' -n' mod N)f {^1,1')- (13) 

l,l' = -N/2 

and the second Riesz transform /i(x„ „/) is obtained analogously by circularly convolving /(x„ „/) with f(^^(x„ „/). 
In practice, if implemented numerically, we always calculate the periodic discrete Riesz transform. This has a direct 
implication on the measure of unidirectionality we propose because it introduces an additional error due to periodic 
filtering. 

III. Second-order statistical characterization of the monogenic signal 

Consider a zero-mean random field /(x) with covariance function r//(x, ^) — cov{/(x),/(x — ^)}, where x 
is a global spatial position and ^ a local spatial offset. If r//(x, ^) — ?'//(^), Vx, i.e., the covariance function 
only depends on the spatial offset ^, then /(x) is called (wide-sense) stationary. Nonstationary random fields are 
difficult to analyze, and most of the properties we will derive do not hold for this general case. Moreover, most 
nonstationary random fields of interest can be locally approximated as stationary for a small enough patch p9) . 
Therefore, from now on we only consider stationary random fields, which is also the kind of random field that 
most of the statistical literature focuses on | [30| . 

In order to characterize the second-order moments of the random monogenic signal, we obtain the covariances of 



the random field and its Riesz transforms. To do so, let us consider the spectral representation |31 1 of the random 
field /(x) |20) 

/(x)=|ydZ;(k)e*'^"^ (14) 

April 11, 2013 DRAFT 



Since /(x) is stationary, the spectral process -Zy(k) has uncorrelated increments, i.e., cov{dZf(k),dZj:{'k')} = 
S'//(k)5(k — k')dk(ik', where SffCk) is the power spectral density of /(x). The covariance function and the 
power spectral density are 2D Fourier transform pairs: 

r-//(0 = //5//(k)e*'^"«dk. (15) 

Taking into account the spectral representation of the Riesz transforms, we may express their covariances as 

Tggi^) = JJ COS^{^)Sff{k)^''"^ dk, (16) 

Thhi^) - J J sm^{K)Sff(k)e'''"^ dk, (17) 

and we note that 

rffiC=rgg{^)+rhh{^). (18) 

This is an important property because it is a relationship between the three covariances, which we can use to 
simplify the statistical description of the monogenic signal. It can be shown that ( 18 i only holds for stationary but 
not for nonstationary random fields. Additionally, we find that the cross-covariances are given by 

rfg{$) = -i J J cos(K)^//(k)e'''"« dk, (19) 

rfh{$) = -i J J sin(^)5/^(k)e"'"« dk, (20) 

rgkiO - I |y'sin(2^)5//(k)e»''"« dk, (21) 

which satisfy 

rgf{i) = -rf9{^). r,,f{^) = -rf^{^), r^gi^) = Vg^i^)- (22) 

This provides a complete second-order statistical characterization of the monogenic signal with six auto- and cross- 
covariances. 

So far we have characterized the monogenic signal in terms of its (real) components, but now we shall obtain the 



quaternion- valued characterization. Section \ll-A\ has shown that in general a covariance and three complementary 
covariances, all quaternion-valued, are necessary to completely characterize a random quaternion. However, since 
the monogenic signal does not have a /c-part, it is obvious that two complementary covariances suffice. As we will 



see shortly, for stationary random fields we only need one complementary covariance due to ( 18 i. 
Let us start with the covariance of the monogenic signal, given by 

TmmiO = cov{m(x), m(x - I)} = 2rff{^) - 2irfg{^) - 2jrfi,,{^). (23) 

This covariance function only specifies three of the six real covariances. In order to access the remaining three real 
covariances, one may be tempted to use the standard complementary covariance. However, it is straightforward 
to show that r,nm*{^) = cov{r7i(x), m*(x — ^)} ~ for all stationary random fields, so this complementary 
covariance does not provide any useful information. Consider instead the covariance between m(x) and m(*)(x) 

rrnmi') (^) = '^rggi^) - 2*r-/g(C) + 2fcrg„(|), (24) 
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(a) Isotropic 



(b) Anisotropic 



Fig. 1: Sample images of random fields 



(c) Unidirectional 



which specifies a further two real covariances. Recalling now that rhh{^) = ^ffiO ^ ^ss(^) gives us access to 
final remaining real covariance. Hence, r„„„(^) and r^^ii){^) together contain the same information as the six 
real covariance functions. We remark that instead of r„j„j(i)(^), one can also use 122J 

r-„^™(.)(0 =cov{TO(x),m(")(x-0}, (25) 

where ?/ is any pure unit quaternion with zero fc-part. The extra freedom of being able to choose an arbitrary basis 
provides flexibility that we shall use later. 

IV. ISOTROPY AND DIRECTIONALITY 

A fundamental characteristic of a random field is the degree of rotational invariance of its second-order statistical 
properties |J9] p. 57]. A spatially isotropic random field exhibits perfect rotational invariance as its covariance 
displays circular contour levels. Fields that are not isotropic are called anisotropic. A common subclass is the class 
of geometrically anisotropic fields [9 p. 61] that have covariances with elliptical rather than circular contours. 
Taking this idea to the extreme, we anive at unidirectional random fields, where there exists a rotation so that there 
is variation only in one of the two axes. Figure [T] shows samples of three random images: isotropic, geometrically 
anisotropic, and unidirectional. In this section we provide formal definitions of these three types of random fields. 

A. Isotropy 

Definition 4.1: A second-order stationary random field /(x) is isotropic if the covariance of the field is finite 



and only depends on the magnitude of the lag, that is if rff{^) = Cj f v^^€) 1 30 . 

Figure[Ta]shows a sample of an isotropic random field. We observe that it is not showing any preferential treatment 
of any given image orientation. That is, if we rotate an isotropic random field, it retains its statistical properties for 
arbitrary rotation angle. We now derive the covariances of the monogenic signal for an isotropic random field. Let 
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the spatial lag be expressed as ^ = ^[cos(0),sin(0)]-^, and let Sff(k) ~ Sff{k) be the power spectral density of 
/(x). Hence, rff{^) is given by 

rffi^) - J J Sff(k)e^^"^ dk = J^ J Sff{k)e^'^"^ kdkdn 

nOO 

= 27T Sff{k)Jo{kOkdk = 27rHo\Sff,0, (26) 

where we have used polar coordinates to solve the integral, Jq() denotes the Bessel function of the first kind and 
order a, and 'H^^(S'//, •) is the inverse Hankel transform of Sff and order a | |32| . Following similar derivations, 
the covariances of the Riesz transforms become 



rggii) = 27T cos^ {e)noHSff,o '^n^\k-'Sff,o, (2?) 

r;,,(0 = 27Tsm'{e)no\Sff,0 + ^-^^nY\k-^Sff,0, (28) 
and, as expected, they fulfill r//(^) = ?'gg(^) + rhhi^)- Finally, the cross-covariances are given by 

rfg{0 = -2ncos{e)n^\Sff,0, (29) 

rf^,i^) - -2iTsmi9)HY'iSff,C), (30) 
rghiO=smi20) 



7:n^\Sff,o--^n^\k-'Sff,o 



(31) 



From the previous equations, it is clear that the six real covariances are functions of only three quantities, namely 
'HQ^{Sff,-), ^{^^{Sff,-) and 'Hi^{k^^Sff,-). Therefore, we may characterize the second-order statistics of a 
monogenic signal for an isotropic random field using one properly selected quaternion-valued covariance, which 
gives us access to the three unknown quantities. For instance, we can use r„„j(i)(^). 

B. Geometric anisotropy 

Definition 4.2: A second-order stationary random field /(x) is geometrically anisotropic if the covariance of the 
field is finite and only depends on the magnitude of the deformed lag as r//(^) — Ca(vs^^)' for some 2x2 
symmetric positive definite matrix D. 



Figure lb shows a sample of an anisotropic random field. We can see that there is a preferred direction, which is 
clearly visible as we have chosen a matrix D with large condition number Now we will investigate the statistical 
properties of the monogenic signal for a geometrically anisotropic random field. The covariance of /(x) is given 
by 

rff{S,)^2^n^\Sff,l), (32) 



where ^ = d (7iS,\ + Cj^^^l, with ai being the largest eigenvalue of D. The proof can be found in Appendix 
Contrary to isotropic random fields, we generally need two quaternion-valued covariances to completely characterize 
the statistics of a geometrically anisotropic random field. To see this, let us consider the covariance of the first 
Riesz transform, given by 

oo 

Z — — OO 
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Fig. 2: Plot of /3(cri, a) for cti = 0.5, and a = and a = tt/2. 



where ai{ai,a) are the Fourier coefficients of 



cos(a) cos(k)\/(7? cos2(k) + sin^f/i) =F sin(a) sin(K)i/cos2(K) + cr, ^ sin^f/t) 

/^('^i. «) = ^ ,, A^ -1 • 2, , ' (34) 

(Ti cos^(kj + (Tj sm (k) 

a is the angle of the dominant eigenvector of D, and the sign depends on the determinant of the eigenvectors matrix 

of D, i.e., whether it is a rotation or a reflection matrix. The proof is also presented in Appendix [A| Hence, just 

considering rgg{^), we see that at least two quaternion-valued covariance functions are required. One could argue 

that some (or almost all) of the coefficients ai{ai,a) might be zero and, therefore, the correlation would depend 

only on a few 'Hi{Sff,£,), as in the isotropic case. To show that this is not true. Figure b^ depicts the function 

f3{(Ji, a) for fji = 0.5, and a = and a = 7r/2. We can see that several Fourier coefficients are needed to express 

P{(Ji,a) even for this simple case. For arbitrary values of cri and a, simulations have shown that there does not 

exist a simple Fourier expansion of the function. 



C. Unidirectionality 

As the condition number of D in the geometric anisotropic model increases, the random field will appear more 
and more directional. However, no matter what finite condition number is selected, the random field is still an 
aggregation of more than one wavenumber coefficient. We therefore introduce an idealized third class of random 
field, namely purely directional signals. 

Definition 4.3: A second-order stationary random field /(x) is purely unidirectional if its covariance takes the 
form rff{^) — Cuin.'^^), where n ~ [ni,n2]'^ = [cos(i^),sin(i^)]-^ is a unit-norm vector in the direction of the 
random field. 
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A random realization of a unidirectional random field can be found in Figure Ic Using the spectral representation, 
a unidirectional random field may be expressed as 

/(x) = ydZ/(fc)e''=""^ (35) 

where cov{dZf{k), dZf{k')} = Sff{k)d{k — k') dk dk' . The unidirectional signal is thus a one-dimensional signal 
embedded in two dimensions |33|. The spectral representation allows us to write the Riesz transforms as 

g(x) = cos(j^)s(x), /i(x) = sin(i^)s(x), (36) 

where 

sgn(fc)dZy.(fc)e*'="'^'' (37) 

is a partial Hilbert transform in direction n |13,|. It is easy to show that the covariances are 

r,g(^) = cos2(i.)r/;(^), rnh{i) = sm\v)rff{i). (38) 

The cross-covariances are 

Tfhii) = i'An[u)rfg[i), rgn{S,) = \ sm{2u)rff{i), (39) 



where 



(^) = / 5//(fc)sgn(fc)e"=""«dfc. (40) 



Similar to the case of isotropic random fields, we do not need five real covariances to characterize a unidirec- 
tional random fields but only two. Thus, only one quaternion- valued covariance suffices, for instance r„„„(^). A 
consequence of this correlation structure is that unidirectional random fields are C-proper: 

Theorem 1 (Characterization of a unidirectional random field): The monogenic signal of a random field /(x) 
is C'-proper, where -q = co&{v)i + sin(i/)j, if and only if /(x) is stationary and unidirectional with direction 
n= [71.1,712]"^ = [cos{v),sv[\{v)Y' . 

Proof: See Appendix [B] ■ 

This result shows an analogy between the monogenic signal of a stationary unidirectional random field and the 
analytic signal of a stationary process, which is complex proper p5) . That is, the complementary covariance of the 
analytic signal x+{t) for a stationary complex process x{t) is zero, i.e., cov{x+(t),2;!j_(t — r)} ~ 0,Vr. 

V. Test for unidirectionality 

In this section, we propose a measure for the degree of unidirectionality in a random field and estimate its 
preferred direction. We then use this measure to build a test for whether there is sufficient statistical evidence to 
classify a random field as unidirectional. 
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A. Measure of unidirectionality 

The monogenic signal of a unidirectional random field may be expressed as 

'7^(x) = /(x) +77s(x), 



(41) 



where s(x) is the partial Hilbert transform in the direction -q of the field. This means that m('')(x) = m{x) for the 
unit quaternion 77 — cos{v)i + sin(z^)j. Therefore, the quantity 

'1 



min£' 



to(x) — ?7l'''-'(x) 



,(0)-maxRe(r„„(,)(0)) 



(42) 



is zero if /(x) is unidirectional and greater than zero if it is not. Based on this mean-squared error, we propose to 

use 

2maxRe(r,„„(„)(0)) 

U = — — 1 (43) 



^7nm [^ ) 

as a measure of unidirectionality. We will see momentarily that this measure is normalized to values between 
and 1. The maximum Z^ = 1 is attained if the field is unidirectional, i.e., m^'^'>(x) — m(x). We now simplify the 
expression ( |43] l for U. We start by considering the numerator. Expressing the involution over 77 in terms of the 
involutions over the canonical basis {i,j, k}, after some tedious algebra the numerator becomes 



maxRe (r„„„(,,) (0)) = max n Rn, 



(44) 



with R 



(45) 



Re (r„„(.) (0)) Imfe (r„„(.) (0)) 
_-Imfe (r„„„(,) (0)) Re (r,„,„(,) (0)) _ 
where Im^ (q) is the fc-component of q. Hence, it is easy to show that ( |43] l can be written as 

, , _ 2Amax(R-) , 

' mm \^ J 

where Amax(R) is the largest eigenvalue R and the direction n is given by the dominant eigenvector. Taking into 
account that the matrix R is given by 

[ 2rgh{0) rff{0)-rgg{0)+rh^,{0) 

we use the closed-form expression for the largest eigenvalue of R to write 



(46) 



U 



is(0) + ^L(O) - 2rgg(0)rw,(0) + 4r2,^(0) 



rff{0) + rggiO)+n,h{0) 
Using the properties of the covariances, U finally becomes 



U 



(47) 



(48) 



2rffiO) 

It is now obvious that U is lower-bounded by 0, so U E [0,1]. To shed some light on U, let us detemiine which 
signals minimize and maximize U. 
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To achieve U = 0, we need ^^^(0) — r^;i(0) = and rgh{0) = 0. This is equivalent to 

ff cos{2K)Sff{k) dk = 0, ff sm{2K)Sff{k) dk = 0. (49) 

Combining both conditions we have 

ff Sff{k)e-'^''dk = 0, (50) 

which may be rewritten as 

/■oo 

/ a2{k)kdk^0, (51) 

where 

02 (fc) = / S'//(k)e-*2«dK (52) 

J —TV 

is the second Fourier coefficient. We would hope that an isotropic random field achieve U — 0. For an isotropic 
random field, the power spectral density is 5//(k) — Sff{k), yielding 

Sff{k)e-'^''dk= Sff{k)kdk e-'^" dn. (53) 

It is clear that the integral in k is zero, and provided that Sff{k) satisfies 

Sff{k)kdk <oo, (54) 



the measure of unidirectionality for an isotropic random field is indeed zero. It may be tempting to think that 
isotropic random fields are the only fields that attain the lower bound of 0. However, from the form of e^*'^" we 
see that to achieve U = a sufficient condition is 

Sffi--ki,k2)^Sff{ki,k2) (55) 

and 

Sffiki,-k2)^Sff{ki,k2). (56) 

Thus, any function that exhibits parity invariance independently in either of the two arguments, e.g., a separable 
function, also leads to W = 0. One may thus argue that, in some sense, separable covariances are as far from being 
unidirectional as isotropic covariances. 

While 14 = is only a necessary but not sufficient condition for a random field to be isotropic, Z^ = 1 is indeed 
necessary and sufficient for a field to be unidirectional. The proof of this statement follows along the lines of 
Appendix IB] 

B. Determining the threshold 

In order to test whether there is statistical evidence to classify a random field as unidirectional, we need to 
determine a threshold for U. This threshold needs to take into account that, in practice, we only observe a single 
finite patch of a random field and therefore work with the periodic discrete Riesz transform. Recall that m(x) is 
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the periodic monogenic signal. We may estimate our measure of unidirectionality from a given realization of the 
random field of size N x N as 

^^2AmaxW_,^ (57) 

^mm [^ } 



where 

Re (ffhM^) (0)) Imfc (f,y„^(.) (0)) 



R 

and 



-Imfc (ffj^fr^u) (0)) Re (f^^o) (0)) 



(58) 



'i^^>W ^ j;p^^(^n,n')m'^^^ (x„,„')> (59) 



and m^^\x) stands for either m(x), m*^*^(x) or m^-'^(x). Then, the estimated measure of unidirectionality may be 
expressed as 

and the estimated preferred direction is given by the dominant eigenvector of R. In order to determine a threshold 
on U, we need to determine the probability of false alarm, which requires knowing the distribution of the statistic 
U under the null hypothesis "the random field is unidirectional". As one may expect, deriving the distribution of 
U is far from trivial. However, for a bandpass random field, we are able to establish the following approximation 
of the false alarm probability. 

Theorem 2: Consider a finite realization of a random field of size N x N whose power spectral density is 
bandpass with lower cutoff frequency A;. Using the estimated measure of directionality U to test whether the field 
is unidirectional, the threshold 1 — 77 may be determined such that the false alarm probability satisfies 



p{u<i-v)<^. 



Nt] 



4 14 14 



(61) 



Moreover, for unidirectional random fields, U ^ 1 as N ^ 00, and the error 1 —U that is due to considering a 
finite patch of the random field decays with 1/A^. 

Establishing these results is very involved and requires lengthy derivations. Readers that are not interested in the 
technical details may therefore wish to skip the following section and continue their reading with Section |VI| 

C. Derivations 

In order to establish Theorem|2J we proceed as follows. We first consider a deterministic plane wave and evaluate 
the error that is due to the finite size of the random field. We then generalize this result to a random plane wave, 
and finally to a unidirectional random field, which we write as an infinite sum of random plane waves. We begin 
with the following lemma. 

Lemma 1: Considering a deterministic plane wave, given by 

/(x) = Acos(A:on'^x + </.), (62) 
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where A, (f>, and fcg — 2ttXq are given numbers, we have 



N/2-1 N/2-1 

U2^l-U^ Y, Y. C{Xi,X2) + 0{l/N^), (63) 

\i = -N/2\2 = -N/2 



where 



C{Xi,X2)-^, 



{Xin2 ~ X2nif 



xi + xi 



sin^(7rA^(Aoni - Ai)) sin^(7r7V(Aon2 - X2)) 
sin^(7r(Aoni - Ai)) sin^(7r(Aon2 - A2)) 



(64) 



with ki — 2TiXi. 

Proof: Similar to f,^j^j(o)(0), we define fff{Q) &P\ 

^//W = ]^E/'(^»^"')- (65) 

The covariance of the periodic discrete Riesz transform g(x) is estimated as 

^5S(0) = J^ E5'(X",«') = rgM + [fgg{G) - rggm = fgg{0) + ^Tgg, (66) 

with 

n,n' n,n' 

and the covariances ff^f^{0) and f^j^{0) are obtained analogously. Performing a Taylor series expansion of U in 

terms of ^rgg, Ar^,~ and ^r^j^, we find that 

2ri9Afss + 2n?AfrT; — 4nin2Af-? / o \ 

U = l- ^ ^^J^ '-^-^ + O (max(Ar ^^ ,^j) , (68) 

r//(0) V (•)(•) y 

where O ( max(Arj^~ - ) J stands for higher-order terms of Ar^g, Ar^^ and Ar~^. We now need to quantify the 
deviation error U2 = ^ —U. Let us start by writing the Fourier transforms of ^(x) and g(x), which are given by 

G(k) = i?(i)(k)F(k) = G+(k) + G-(k) (69) 

and 

G(k) = i?(i)(k)F(k) = G+(k) + G-(k), (70) 

respectively. Here, G* (k) and G* (k) are the infinite-length discrete-space and periodic discrete first Riesz trans- 
forms corresponding to the positive and negative frequencies, respectively. These are given by 

G^fk) = in ^c^'^ ^ («-i)[(fco"i=Ffci)+(fco^2Tfc2)i sin(iV/2(fconi =F fci)) sin(iV/2(fcon2 T h)) .^^^ 

^ '^ ^'^i 2 e e sm{l/2{koni T h)) sin(l/2(/cori2 T ^2)) ' 

^±..^_ . fci A ^^ ^ (iv-i)[('.o^iT'=i)+(fco"2=Ffc2)i sin(A^/2(fconi =Ffci))sin(A''/2(A;on2Tfc2)) „,. 

^^)- '^^^rq2^ "" ' sin(l/2(fconiTfci)) sin(l/2(fcon2TA;2))' ^ ^ 

We notice that the main difference between G*(k) and G*(k) is the first term, which is ni for the former and 
ki/y'kl + fc| for the latter The term ni comes from obtaining the Riesz transform of a complex exponential, 

^There are similar definitions for rgg(x), r(i(i(x), and fgi^{K.). 



and 
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which may be done in closed form. The term k^j \Jk\ + k\ is precisely the filter that has to be applied to obtain a 
Riesz transform. Now, the Parseval-Rayleigh relationship allows us to express ^^^(O) in the frequency domain as 



ki,k.2 



fci,fc2 



f+M + f^M+Oil/N'), (73) 



and similar expressions hold for fgg(O), '^^^(0) and ?^5^(0). Therefore, the error term for the positive and negative 
frequencies becomes 



99 gg^ ' 99^ ' 4^4 Z^ 

Ai,A2 



A? 



xi + \l 



s\T?{TTN{XQni T Ai)) sin^(7rA''(Aon2 T A2)) 
sin^(7r(Aoni T Ai)) sin^(7r(Aon2 T A2)) 



(74) 



where we have used fc, = 27rA,, and the overall error becomes 



Af^^ = Aft + Af,- = ^74 E 

Ai,A2 



99 2N^ 



A? 



A^ + Ai 



sin^(7r7V(Aorii - Ai)) sin^(7riV(Aon2 - A2)) 
sin^(7r(Aoni - Ai)) sin^(7r(Ao?i2 - A2)) 



0{1/N^). 



(75) 



Applying the same procedure, we obtain similar expressions for the remaining error terms. The error in the covariance 
of /i(x) becomes 



Af 



hh 



2N^ 



E 

ki,k2 



Ai 



Af + Ai 



sin^(7r7V(Aoni - Ai)) siii^(7rA^(Aon2 - A2)) 



and for the cross-covariance it is 

A1A2 



Ar.r, = :r^ E 



gh 2A^4 



Ai ,A2 



A? + Ai 



sin^(7r(Aoni - Ai)) sin^(7r(Aon2 - A2)) 

sm^{TTN{Xani - Ai)) sin^(7r7V(Ao7i2 - A2)) 
sin^(7r(Aoni - Ai)) sin^(7r(Aon2 - A2)) 



+ Oil/N^), (76) 



+ 0{l/N^). (77) 



Finally, the proof is concluded by inserting ( fTS] ), ( |76] l and ^T7\ into ( |68] l and taking into account that r//(0) — A'^/2. 



Let us now analyze ( |63| ) more carefully. To do so, we decompose XqUiN as XqUiN — lXoniN~\ + ci, where 
[XqUiNI and c/ are the integer and fractional parts of XqUiN, respectively. Defining the transformed indices 
ji = ki — {XqUiN^, we may write C(Ai,A2) = C{ji,J2)- We can consider different regions for the orders of 
magnitude of the transformed indices and analyze C(ji,J2)- For small values of both indices, i.e., ji — 0(1), 
the function becomes C{ji,J2) ~ 0{l/N'^). For large values of only one of the indices, i.e., ji = 0(1) and 
ji = 0{N) or vice-versa, we find that C{ji,J2) = 0{1/N'^). Finally, for large values of both indices ji = 0{N), 
we have C(ji,j2) = 0{1/N^). Now, defining the sets J^a' = {{J1J2) ■ ji = 0{N"),J2 = 0{N°'')}, we may 
decompose U2 as 

t/2= 601 +610 -600 + 611, (78) 



where each term is given by 



Baa' = E C'(ji,J2) 



(79) 
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Finally, taking into account that the size of the sets J^a' is \Jaa' \ — A^"+" , it is easy to prove that Sqi = 0(1/A^) 
and Sio = ©(l/iV), whereas 600 = ©(l/iV^) and Su = ©(l/iV^). Hence, the error is 

C/2 = 601 + Sio + 0(l/iV2) = 0(1/7V), (80) 

that is, it decays linearly with N . Hence, from the previous analysis we find a simpler expression for U2, which is 
given in the following lemma. 

Lemma 2: For a deterministic plane wave, the error due to the finite periodic discrete Riesz transform is 

U2 = ^ [sin'(^C2)5(Ao,ni,n2) + sin2(7rci)g(Ao,n2,ni)] +0{l/N^), (81) 

where Q{\a,ni,n2) = g+{\o,ni,n2) + g^{Xo,ni,n2) and 

r.1/2 1 ^2 



f ' 1 A^ 

Jo Agnf + (Aon2 ± A)-^ sm^ (ttA) 



(82) 



The convergence rate with which U ^ I as N —>■ 00 is therefore linear in N. 
Proof: Let us consider Bqi, which is 

2 "^^^ A {{ci/N + j,/N)n2-{c2/N + J2/N)mf 



\n\=^' 31=--^ 



sin (7r(ci + ji)) sin (7r(c2 + J2)) 



-, (83) 



sin^{7r{ci/N + ji/N)) sm^{TTic2/N + J2/N))' 
where L ^ N and /3 is a real number close to one, but smaller Taking into account ci,ji ^ N and sin(x) « a: 
for small x, we find that 

J_ ^ sin^(7r(ci+ji)) ^^' (c2/iV + j2/JV)^n^ sin^(^(c2 + J2)) 

It is easy to verify that for a relatively small value of L, the first sum in the above expression is almost one for 
every value of ci. Moreover, sin^(7r(c2 +.72)) = sin^(7rc2), which yields 

B ^ A "v^' i^2/N + j2/NYnl sin^(7rc2) 

°' ^ ^' |,,Pv^ ^0"? + (^o"2 + C2/iV + j2/iV)2 sin2(^(c2/7V + j2/A^)) ■ ^ ' 

The above sum may be considered as a Riemann integral, and therefore 

Bo, « ^nlsin\nc2) J^^^^^^^_^ A^n? + !aoC+ c'/iV + A2)^ sm^(.(c2)A^ + A2)) ''^- ^''^ 

Applying the change of variable A = C2/N + A2, and observing that N^^^ w 0, we find that 



2 
■Boi « — nisin^(7rc2) 



1/2 ^ ;^2 .1/2 ^ ;^2 ■ 

«A + / . „ „ 7- --TT 7T-, "A 



Agnf + (Aon2 + A)2 sin^(7rA) Jq A§n^ + (Aon2 - A)^ sin^(7rA) 



(87) 

Finally, the proof follows from the symmetry of Bio with respect to Boi- ■ 

Before considering random plane waves, we will obtain an approximate closed-form solution to the integral in 
(|82jl, given in the following lemma. 
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Lemma 3: The integral Q±{Xo,ni,n2) is approximately given by 

^2 / 1 A, 



g±{Xa,ni,n2) ~ y 



I'^il 



^0 



Ao^2 



{2nj - I] 



arctan 



1 , n2_ 
2Ao|ni| \ni\ 



=p arctan 



n2_ 
\ni\ 



T 



Aonfn2 



"■2 



1 



"•'^'^t„Um 



(88) 



Proof: As far as we know, there is no analytic solution for the integral ( |82| l. Hence, we may consider the 
second-order Taylor series only of the second term of the integrand, and thus 

1 / TT^X^' 



t/±(Ao,ni,7i2) 



1/2 



which may be split as 



g±{Xo,ni,n2) 






1/2 



A2 ± 2Aon2A + A§ 



dX 



1 



1/2 



dX, 



A2 



dX. 



(89) 



(90) 



A2±2Aon2A + A^ 3 J^ A2 ± 2Aon2A + A§ 

By applying a partial fraction expansion to the integrand of the second integral, we may find the primitive functions 



e.g., in 1 34 1, and by substituting the integration limits the proof follows. 

In U2 we need the sum of 5+(Ao,ni,n2) and t^_(Ao, ni, 712), which may be simplified as follows. 
Lemma 4: tJ(Ao, ?^i, ?^2) is given by 



^(Ao,ni,n2) = (y+(Ao, 77,1,77,2) + 5_(Ao, 771,772) 



7^1 Ao77i772 



loe 



/ 4Ao772 \ 

i+4Ar 

4Ao772 

v^ + TT^y 



|77,l 



1 



Proof: It is easy to show that 

t;(Ao, 771,712) = C;+(Ao,77l,772)+t/_(Ao,77l,772) 
1 Ao 



Ao^2 
77i XQn\n2 



y (2«2 



arctan 



/ 4Ao|77i| 

Vxi-\ 



(91) 



I'^il 



277^ - 1 



3 

arctan 



Now, using 



,2Ao|77i 
arctan (a) + arctan (6) = arctan 



log(^l 

_ »2 



772 

Ao 



4A2 



arctan 



(mod 7r) 



- log ( 1 
1 



772 

Ao 



2Aobi| 



772_ 
77l| 



4A2 



(92) 



(93) 



\-ah . 

the proof follows. ■ 

So far, we have only considered a deterministic plane wave. In the following, we will generalize our results to a 
random plane wave with random amplitude, phase, and direction but fixed frequency. If we would like to test the 
null hypothesis that a realization of a random field is a plane wave vs. the alternative that it is not unidirectional, 
we need to choose a threshold for Z^ (or C/2 = 1 — W). In order to determine the probability of false alarm, we 
would need to know the distribution of U2 under the null hypothesis. As deriving this distribution seems difficult, 
we instead obtain a conservative bound on the probability of false alarm using Markov's inequality p5] . Concretely, 
Markov's inequality states that 



P (C/2 > 7?) < 



m2 



(94) 
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for a positive random variable U2- Hence, we have to find £[1/2] for a random plane wave, which is presented in 
the following theorem. 

Theorem 3: The expectation of U2 for a random plane wave is 

1/41 4, 1 4 \ 

Proof: The proof is presented in Appendix [Cj 
Now, using Theorem [3] and Markov's inequality, we may bound the probabiUty of large U2 as 



(95) 



1 



^(^^^^)^]v;^l^A^-^° + ^ 



4 1 



(96) 



which is only valid when the right-hand side is smaller than one. The following theorem generalizes this result to 
unidirectional random fields. 

Theorem 4: The expectation of U2 for a unidirectional random field is 



E[U2 



1 



1/2 



4 14 1 



S{X) dX 



N 



1/2 



Oil/N^). 



(97) 



S{X) dX 



Proof: The proof is presented in Appendix [D] ■ 

This shows that, for unidirectional random fields, the rate with which _E[Z^] ^ 1 as A^ ^ 00 is 1/A^, just like 
it was for deterministic plane waves. Unfortunately, E[U] still depends on the power spectral density (PSD) of the 
random field. Since this quantity is generally unknown, we should bound E[U] to obtain an expression that does 
not depend on the PSD. In the derivations for the random plane wave, we already required that the frequency of 
the plane wave was not too low or too high. We shall therefore assume that the random field is bandpass. 
Lemma 5: The expectation E[U2\ may be bounded as 



EP2] < 



4 14 1 



-Oil/N'), 



where A; is the lowest frequency. Equality in ( [98| l is attained if the PSD takes the form 

S{X) = <j^S{X-Xi). 



(98) 



(99) 



On the other hand, if the PSD has support at higher frequencies, then ( |98| l is a strict inequality. 
Proof: Let us define 



S{X) 



5(A) 



(100) 



S{X) dX 



where A/ and A^ denote the lowest and highest frequency component of the PSD. Observing that S'(A) integrates 
to one, ii^[C/2] may be written as^ 



E[U.]^- 



4 1 4 \ ~ 1 



(101) 



*For the sake of notational simplicity we ignore tlie term 0{1/N^). 
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Now consider the integral 



Ai 



4 1 4. 



X]Si\)d\. 



(102) 



By taking into account that the function in parentheses is monotonically decreasing and positive on the considered 
interval, we may write 



^^ '^x)s{X)dX<(A,^-^\i] I 5(A) dA 



4 1 



A/ 



7r2 A, 9 " 



(103) 



IX, \TT^ X 9 J ' ' - \tt^ Xi 9 
and the proof follows. 

Theorem l2] is now obtained by applying Markov's inequality 



VI. Numerical results 

In this section, we present simulation results illustrating the behavior of our measure of directionality, which we 
analyzed theoretically in the previous section. First, we examine the statistical behavior of our measure for different 
kinds of random fields. Then, we apply our measure to a real-world problem, where we detect unidirectional patches 
on the surface of Venus. 



A. Statistical behavior of LI 

We generated samples of isotropic, geometrically anisotropic, and purely unidirectional random fields. The PSD 
of the isotropic and geometrically anisotropic random fields is a shifted (band-pass) Matern covariance function 
pO| , given by 

SffW = 



a^T{v+l){4iy)'' 



7rr(t/)(7rp)2'' 



4j./(7rp)2+(VA^"DA-Ao)^ 



v+l ' 



where k = 2tt\, r() denotes the Gamma- function, and D = I2 for the isotropic field, and 



D 



1 



1 0.85 
0.85 1 



i^+i 



V0.2775 
for the anisotropic random field. The purely unidirectional random field has PSD 

^ff[^) ~ r : 

nr{i^){npy- 4i//(7rp)2 + (A - Ao)' 

along the line A = An. For all fields, we chose ly = 1.5, p = 20, a^ = 1 and Ao = 0.1 as the center frequency. 

We computed our estimated measure of unidirectionality U and estimated the probability density functions (PDFs) 
of U for the three different kinds of random fields of size 64 x 64 from 10, 000 reaUzations. The results are shown 
in Figure pi We can see that the isotropic random fields have U close to zero, and the unidirectional random fields 
have U close to one, in some cases even identically one. These cases happen when the direction is aligned with 
one of the coordinate axes, i.e., u = 01 i^ = 7r/2. The PDFs of isotropic and unidirectional random fields are 
clearly separated and, while the PDFs of anisotropic and unidirectional random fields are closer, they are still easily 
distinguished. For larger values of N, the PDFs are even more clearly separated. 



April 11, 2013 



DRAFT 



20 



12 
10 



to 

o 



■ Unidirectional 

■ Isotropic 

■ Anisotropic 




Fig. 3: Estimated PDF of the three different kinds of random fields with N ~ 64. 



■ Monte Carlo 

■ Theoretical approximation 




Ao 



Fig. 4: Expectation of U for a random plane wave. 



Next, we would like to corroborate the results of Section V-B Consider first a random plane wave, with fixed 
amplitude, random phase and direction. Figure UJ shows the value of Epl] predicted by Theorem pi and the estimate 
obtained from Monte Carlo simulations, as a function of the frequency Aq. We can see that there is good agreement 
except at low frequencies. When Ao becomes small, e.g. when Aq — 0{l/N) the higher order error terms in 
Section V-C become appreciable, see e.g. Lemma [3] The intuitive understanding of this is that when the image 
contains variation over longer spatial scales than the patch we have observed, this will inevitably produce problems 
with our estimation methods. We therefore assume we have observed a large enough patch to characterize the 
important structure of the random field. This is a very intuitive effect. 

Finally, we consider unidirectional random fields, which are generated as explained at the beginning of this section. 
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N 

Fig. 5: Expectation of 1 —ZY for a unidirectional random field. 



Figure [sj shows the expectation of 1 — Z^, obtained through Monte Carlo simulations, the theoretical approximation 
given by Theorem |4] and the upper bound given by ( |98| ). For the upper bound, we have assumed that the lowest 
frequency component in the PSD is at A; = 0.05. We can see that the theoretical approximation is very good. While 
the bound is not very tight, it does display the same behavior with respect to the size N. That is, it decays with 
l/N, just like the Monte Carlo simulations. By using an estimated PSD in Theorem [4] we obtain a tighter bound, 
which is shown in magenta. 



B. Application 

We now present an application of our measure where we would like to detect unidirectional patches on the surface 
of the planet Venus. The image we consider is depicted in Figure [6] We would not expect the entire image to be 
unidirectional, but there are clearly some unidirectional patches. We therefore apply our measure of unidirectionality 
in a sliding-window fashion. For each pixel in the image, we calculate U for a 16 x 16 neighborhood centered around 
that pixel. Three such neighborhoods are marked with a square in Figure l6] and shown in greater magnification in 
Figure It] The arrows in Figure [6] indicate the estimated direction of the three patches. 

Looking at the three patches in Figure [7] in detail, we see that Patch 1 (which corresponds to the bottom left 
square in Figure l6]), is almost unidirectional and Patch 3 (which corresponds to the bottom right square), has a 
strong unidirectional component. Patch 2, on the other hand, is the least unidirectional of the three patches. These 
observations match the estimated degrees of unidirectionality. We see that the automated measure can therefore be 
used to process large volumes of data automatically - clearly recognized regions of anisotropy can then later be 
put to manual scrutiny. This is of strong interest in Earth Science, see e.g. 136). 
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Fig. 6: Topography of Venus. 



VII. Conclusions 

In this paper, we have introduced the random monogenic signal as a means to analyze directional structures in 
images. Our main contribution compared to previous approaches is the definition of a measure of unidirectionality, 
which may be used to test an image for the presence of unidirectional components. There are two main possible 
extensions of our work. First, we were only able to derive a rather loose bound on the detection threshold, so a 
tighter bound would be desirable. Second, directional structures are usually limited to smaller patches, and random 
fields may exhibit directionality only in certain frequency ranges in the wavenumber domain. Thus, a space- and 
wavenumber-localized version of our test should be the focus of future research. 
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(a) Patch 1: W = 0.906 



(b) Patch 2: W = 0.412 



(c) Patch 3: W = 0.784 



Fig. 7: Several patches of Venus, as outlined in Figure 6 



Appendix A 
Derivation of the covariances of a geometrically anisotropic random field 



The covariance function of a geometrically anisotropic random field is r-ff{^) = Ca(\/^"^D^), so its power 
spectral density may be expressed as 5'//(k) = Sjf ( •\/k^D~ikj. Let us find an alternative expression for 

.„(0 = //5„(v^^D=Tt)..^'..k. (104, 

Using a square root matrix of D, which satisfies D = D^/^D^/^, we introduce the change of variables 

k' = D"i/^k, =» k = D^/^k', ^dk= \det(T>^/^)\dk' ^ dk', (105) 

where D^^^^ = S^'^/^U-^, with S and U the matrix of eigenvalues and eigenvectors of D, respectively. Therefore, 
the covariance becomes (with some abuse of notation) 



rffi^)- SffiWme''' ""^dh'^ Sff mi) e^^'^"'^dk 



Letting ^ = B'^/^^ = US^/^^ = [|cos(^),|sin(^)]^, it is possible to write r//(|) as 

ry^(^) = / / Sff (fc) e'''^''°'^^'^^kdkdK = 2TTHo\Sff,i). 
Jo Jo 

The covariance of ^(x) is given by 

rggii) = JJcos{^)Sff (Vk^D-ik) e'l^'^rfk, 
and considering the change of variables 

k' = U^k,=» k = Uk',^ dk = |det(U)|dk' = dk', 

we find that 

rggd) = f f cosiK)Sff (Vk'^S-ik') e'^'^^^^dk'. 



(106) 



(107) 



(108) 



(109) 



(110) 
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We would like to eliminate cos(k). To this end, we need to rewrite cos(k) as a function of the new variable. Any 
2x2 orthogonal matrix may be expressed as either 



Ui 



cos(a) — sin(Q;) 
sin(Q;) cos(q!) 



or Uq 



cos(q:) sin(a) 
sin(Q;) — cos(a) 



where Ui represents a rotation and U2 a reflection. Thus, we find that 



cos(k) 
sin(K) 

cos(k) 
sin(K) 
The covariance then becomes 



= Ui 



= U2 



cos(k') 
sin(K') 

cos(k') 
sin(K') 



cos(q!) cos(k') — sin(Q!) sin(K') 
sin(Q;) cos(k') + cos(q!) sin(K') 

cos(q!) cos(k') + sin(Q!) sin(K') 
sin(Q;) cos(k') — cos(q!) sin(K') 



(111) 



(112) 



(113) 



(I) = cos(a) / / cos{K)Sff (Vk^S-ik) e'^^^^^dk 



' gg\ 

T sin(a) If siii(K)5'// (Vk^S-ifc) e'^^^^^dk, (1 14) 

where the negative sign corresponds to the rotation matrix and the positive sign to the reflection matrix. Now, let 
us apply the change of variables 



ik'^Si/^U^ 



«rfk' 



k' = S-i/^k, =^ k = S^/^k', ^ dk = |dct(Si/2)|dk' = dk' 
to obtain 

rggiO = cos(a) / / cos(«;)5'// (||k'||)e' 

TMo^)JjM^)Sff{\\k'\\)e^- 

where we notice that this expression also involves cos(k) and sin(K). These are given by 

fci a^ k[ cos(k') 



(115) 



,nn ..k-sV2u-^^j^, 



cos(k) = 
sin(K) = 



^2 (72 ^2 sin(K') 



VW+^ ^aik'^ + a2k'^ y^CT2^os2(K')+sin'(K') 
so the covariance finally becomes 

rgg{^) = jj P{ai.a)Sff (fc) e*"«dk, 

where 



fi{<Ti,a) 



cos(q;) cos{k)J (jI cos2(k) + sii<?(k) =f sin(a) sin(K)Wcos2(K) + a^^ s\v?{k) 



(Ti cos2(k) + ctj^ ^ sm^{K) 

Let us consider the Fourier series of /3(Ai, a), given by 

00 

/ — — 00 



(116) 



(117) 



(118) 



(119) 



(120) 



(121) 
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where the Fourier coefficients are 

aiiai,a) = ^f p{ai,a)e'''^ cIk. (122) 

271- 7_^ 

Note that we cannot guarantee that the number of non-zero Fourier coefficients is finite for any given choice of Ai 

and a. Plugging the Fourier series into the covariance, we get 

rgg{^)= ^ ai{cri,a) Sff{k)kdk e'''^+*'^«™"('=-«)rfK. (123) 

Finally, considering the change of variable 

k' = -+£,- K,^ K^ - +i- k',^ dK^ dn', (124) 

we may rewrite the covariance aaj 

QQ _ POO n-K 

rgg{^)^ ^ i'e'^^ai{ai,a) Sff{k)kdk ^^k^^in(^')-^i>-' dn' 

l=-oo -^0 ■'-'' 

QQ _ /-OO 

^ i'e*%(CTi,a) / Sff{k)Ji{ki)kdk (125) 



/ — — oo 



which yields 



rggiO^ E ^'e^'^Cii{<j,,a)ni{Sff,0- (126) 



Appendix B 
Proof of Theorem[T] 

If /(x) is unidirectional, we have seen that its monogenic signal admits the representation 

m(x) = /(x)+7?s(x), (127) 

where ij — cos(z^)i + sin(t/)j. Define the second quaternion rj' = — sin(i^)i + cos{i')j, and augment the system 
with k. Then 

m(''')(x) = m('')(x) = TO*(x). (128) 

We know that any stationary random field satisfies 

cov{m(x),m*(x-|)} = 0, (129) 

and so the sample from the monogenic unidirectional signal is C -proper 

To prove the converse of the statement we assume that we have the monogenic signal of a stationary random 
field, and that any finite sample from it is C'-proper Because it is a monogenic signal it has the form 

m(x) = /(x) + i.g(x)+jMx), (130) 

^Note that we do not need to change the integral hmits due to the periodicity of the involved functions. 
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without a fc-component. Signals that are C -proper satisfy 

^™r„(V) (^) = cov {to(x), m(''')(x - ^)} = 0, r„„(., (^) - cov {m(x), m('=)(x - ^)} = 0. (131) 

The second statement is true because we have a stationary monogenic signal, and therefore provides no additional 
information as ?7i*-'^^(x) = m*(x). We write the monogenic signal in the basis {ri,r]',k} as 

m(x) = /(x) + 7] {cos{u)g{x) + sin{u)h{x)) + r( (— sin(i^)(7(x) + cos(i^)/i(x)) . (132) 

We therefore find that 

Vm(V)(0 = r^f{i) - cos(2i.)[r,<,(^) - r„;,(|)] - 2 sin(2i.)r,^(|) 

+ r;'[sin(j/)2r/g(^) - 2 cos(i.)r/ft(|)] + /c{[rgg(^) - r,,,,(^)] sin(2i/) - 2rg,,(|) cos(2z/)}. (133) 
Since we are assuming C-propriety, the complementary covariance satisfies ?"„„(»,') (^) = 0, which yields 

tanz.= ':^, (134) 

r,./.(^)=sin2(i.)r//(0, (135) 

rg,(^)-cos2(^)r//(|), (136) 

rgyXi) =sin(i^)cos(j^)r//(^). (137) 

Consider a finite patch of the random field and assemble the samples of the random field and its Riesz transforms 
into the vectors f , g = a^f + hgi^ and h — ayS- + ^hS-t ■ ^hc covariance matrices of the vectors g and h are 

cov{g,g} = a2^ov{f,f} + 62cov{fg^,fg^} =cos2(t.)cov{f,f}, (138) 

cov{h,h} = a^cov{f,f} + 6^cov{f;^,f^} = sin2(i/)cov{f, f} , (139) 

which implies that the covariance matrices satisfy 

cov { f^^ , fg^ } = cov { 4^ , f,^ } = cov {f , f } , (140) 

and the scalars satisfy a^ + 5^ = cos^(i^) and a\ + b^ = sin^(i^). Taking into account the covariance matrices 
between f and g, and between f and h, it is straightforward to show that ag ~ ah ^ 0, which implies that 
6g = cos^(i^), and h\ = sm^{v). Now, using the covariance matrix between g and h, we have that bg = cos(z^), 
and bh = sin(i^). We shall call f^ = s. Then, the Riesz transforms are g = cos(i^)s, and h = sin(i^)s. As we 
can write the sampled monogenic signal in this form for arbitrary sample length, we have g{'x.) — cos(j^)s(x), and 
/i(x) = sin(i^)s(x), or in the frequency domain 

G(k) = -i cos{k)F{\s.) = cos(i^)S'(k), (141) 

i7(k) = -isin(K)F(k) = sin(i/)S'(k). (142) 
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Finally, for this to hold for all k, we must have 

F(k) (xS{k- {iy±rmT)) 

S'(k) = -isgn(n^k)F(k). 
This shows that the only monogenic stationary random field that is C'-proper is the unidirectional field. 



(143) 
(144) 



Appendix C 
Proof of Theorem[3] 



In this appendix, we obtain 



E[U2] = £^[sin^(7rc2)t;(Ao,ni,n2)] + ^[sin^(7rci)t;(Ao, ^2, ni)], 



(145) 



where the expectation is with respect to v, with n = [rii, ^2]^ = [cos i/, sin v]'^. We first note that we may assume 
that Ci is a random variable uniformly distributed on [—1/2, 1/2] and independent of ly, which is approximately 
true for large values of N. Then, £[1/2] becomes 



E[U2] = E[sm\TTC2)]E[g{Xo, ni, 712)] + ^[sin2(^ci)]£;[g(Ao, ^2, 
= -E[g{Xo, Til, ^2) + g{Xo, n2,ni)]. 



nij 



(146) 



Plugging in the values of g{Xo,ni,n2) and Q{XQ,n2,ni), this becomes 

/ _ 4Aon2 \ 
' l + 4Ag 

4Aon2 ' 6 



E[U2 



1 , Ao 



6 6 



E 



1 



n^ 712 log 

Ao 



V 1 + 4A2/ 



Ao 



71,2 ?ll log 



/ _ 4Aoni \ 

' 1 + 4A^ * 

, 4Aoni 



V 



2Ao7r2 6 

^ Ao 

2Ao7r2 6 

1 Aovr 



[rill arctan 



1 71.2 1 arctan 



EWn, 



XqTT 



4AoM 
4A^-1 

4AoH 
4A2-1 

E[\ni\n 



Ao 
3 

Ao 
3 



E 

E 
1 



|rji|n2 arctan 



i + 4A2y 

/4Ao|ni 



V4A^-1 
/4Ao|n2| 



[712 [Til arctan 1 . , 2 

\4Ao - 1 



AoTT 



E[\n2\ 



XqTT 



E[\n2\nl]. (147) 



^2Ao7r 6 J " " 3 '' ' ^' \2Xqtt 6 J " " 3 
To prove the theorem we need to find each of the above expectations w.rt. 1/, which we will take as uniformly 
distributed between and tt. 

Lemma 6: The value of the first expectation is 



E 



n\n2 log 



l + 4Ag 

4Ao772 

V l + 4A2y 



- ^Ao + -Aq. 



(148) 
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E 



Proof: The expectation may be rewritten as 
/ 4Aon2 \ 

n\n2 log 



l+4Ag 
4Aon2 



V 1 + 4A2; 



cos (t^) sin(z^) log 



/ 4Ao sin(;/) \ 
1 + 4A§ 
4Ao sin(;/) 
V + 1 + 4A§ / 



dv 



An 



sin(i/) log 



/ _ 4Aosin(^) \ 
1 + 4A§ 
4Aosin(^) 
V + 1 + 4A2 / 



diy 



An 



sin(3;/) log 



h 



For the sake of notational simplicity let us define a = 4Ao/(l + 4Aq). We may therefore write 



h 



An 



■ f -,, ,'l-asm{v)\ 
sm(i/) log I z : — ^^ I di/, 



and, using for instance |34|, it is easy to show that 

1 



/i = 



2a 



1 + asin(j/) 
Vl - a2 - ij = -Aq. 



On the other hand, we have 



whose solution is 



l2 = 



An 



sm 3t/ log — r—- dv, 

' 1 + asm(i/) 



i(--^-^^7!^)^^«' 



which concludes the proof. 

Lemma 7: The value of the second expectation is 



E 



n^rii log 



l + 4Ag 
, 4Aoni 



-Aq + -Aq. 



V 1 + 4A2/ 

Proof: The proof follows along similar lines as the previous lemma 
Lemma 8: The value of the third expectation is 

/4Ao|ni| 



E 



\ni\ arctan 



\AXl-l 



= -2Ao. 



Proof: The expectation is given by 

E 



I , /4Ao|ni| 

Til I arctan -—^ 

\A\„ ~ 1 



I / M / 4Ao cos(i') 

cos(zy) arctan -^5 

* 4A^ - 1 



dv. 



It is clear that the integrand is symmetric with respect to n/2, which yields 



I / \i /4Ao| cos(t/)| 

I cos(j^)| arctan I — — 2 

\ 4Ao - 1 



dv^'r\osiv).rct.n('^^ 
^Jo ^ ' V 4Ag-l 



/ ^ 4Ao sin(^) N^ 
l + 4Ag 
4Ao sm{v) 
V^ 1 + 4A2 / 



h 



or equivalently 



--»..„e.„(lW^^. ^ 



7r/2 



dv. 



cos(z^) arctan (6 cos(z^)) dv, 



dv . 



(149) 



(150) 



(151) 



(152) 



(153) 



(154) 



(155) 



(156) 



(157) 



(158) 
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where b = 4Ao/(4Aq — 1). Now, using |34|, the integral becomes 



-/2 1 / . X 

cos(i^)arctan(6cos(i^)) dv ~ -- iyb"^ + 1 — 1] 



and substituting the value of 6, the proof follows. 
Lemma 9: The fourth expectation is given by 



E 



\ni\nn arctan -—t. 

* 4A2 - 1 



'^0 2 3 



Proof: The expectation is given by 

E 



ni|n^ arctan ^'^°""^> 



V4A2-1 



1 r, / M ■ 2/ N f4.Xo\cos{v)\\ , 

— / I cos(;/)| sm (z^j arctan -k | dv, 

TT ./n V 4An - 1 



and taking into account the symmetry around tt/2, it may be rewritten as 

/4Ao|ni 



E 



\ni\nn arctan , , ,„ 
' ' ^ V4Ao - 1 



1 





it/2 



2 r^^ / X . 2/ X f4XocosM\ , 

— / cos(i^j sm (i/j arctan „ | di' 

TT Jo V 4A^ - 1 



/4AoCos(j/)\ 1 

cos(i^) arctan — -^5 ] di' 

27r ./o ^ ^ V 4A2 - 1 ; 27r ^0 



"/' ,^ , , /4AocosM 

cos(3i^) arctan ^ 

V 4A2-1 



The solution to Is follows from the previous lemma, and the solution to I4 is 



1 



12&3 



36^ 



562 



which follows from p4| . 
Similarly, we find that 



E 



1 712 1 arctan 



/4AoK[ 
/4AoK[ 
The last expectations are given by 



E 



\n2\n1 arctan 



= E 
^E 



sin(j/)| arctan 



VFTT 



4Ao|sin(;y)| 
4A2-1 



= -2Ao, 



(159) 



(160) 



(161) 



dv . (162) 



sin(j/)| cos2(z/) arctan 



4Ao|sin(t/)| 
4A2-1 



Ao 2Aq 
T^ ~3~' 



E[\n^]=E[\co^{v)\] = -, 

IT 

E[\m\nl] ^E[\cos{v)\s\n\v)\ - |-, 



-E;[|n2|n?] =^[|sin(i/)|cos2(i.)] = 



37r' 



Finally, taking into account all individual integrals and after a lot of tedious algebra, the proof follows. 



(163) 



(164) 
(165) 

(166) 

(167) 
(168) 
(169) 



Appendix D 
Proof of Theorem|4] 

We shall start with the Fourier transform of /(x), which is 



F(A) = 5]/(x„ 



-2iriA-' x„ 



(170) 
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where k — 27rA. Using the spectral representation of the random field, it becomeaj 

F(A) = 



^^^ dZ (A)e"(^-i)[(^"i-^i)+(^"^-^=)l sm{Nn{Xni - Ai)) sin(A^7r(An2 - A2)) 



dZf{X)eJ 
-1/2 

F+{X) + F-{X) 



sin(7r(Ani — Ai)) sin(7r(An2 ~ ^2)) 

n{N-l)[{\n^-\^)+(\n2-\2)] sin(A^7r(Ani - Aj)) sin(iV7r(An2 - A2)) 

sin(7r(Ani — Ai)) sin(7r(An2 — A2)) 



(171) 



Moreover, it will be useful to know (to easily compare with fgg{0) and fhfi{0)) 



E[fffiO)]^rff{0) = ^E 



E/'(^".»') 



7V4Z^| 



E 



|F+(A) 



E 



\F-{X) 



E [i^+*(A)i^-(A)] + E [i^+(A)F-*(A)] | , (172) 



where we have applied the Parseval-Rayleigh relationship, and the expectations, unless otherwise stated, assume a 
fixed direction. Taking into account that the spectral process is proper, i.e.. 



cov{dZf{X),dZf{X)} = Sff{k) dX dX6{X - A), 
cov{dZ/(A),dZ/(-A)} = cov{dZ/(A),dZ;(A)} = 0, 
where we have used the fact that /(x) is real, we find that 



E 



\F-{X)\ 



= E 



\F^W\ 



and 



This yields 



E [i^+(A)F-*(A)] = E [F+*(A)i^-(A)] = 0. 



Eirffm-^T. ^w^^ 



1/2 



^ sin^{TTN{Xni - Ai)) sm^{TTN{Xn2 - A2)) 
sin (7r(Ani — Ai)) sin (7r(An2 — A2)) 



Let us continue with the periodic discrete Riesz transform 



'"^* ^ -'TWrn'"'^ ^ -'7WTM [^'''' " ^"'"'i 



and its square 



G(A) 



A? 



xi + xi 



\F+{X)+F-{X)\\ 



which allows us to write 



'^-(0) = 42E5'(-"-') = ^E|g(A) 



' ggy 



A 

\2 



1 V Af 

m^ xi + XI 



(\F+iX)\'' + F+{X)F-*iX)+F+*{X)F-{X} + \F-{X) 



(173) 
(174) 

(175) 
(176) 
(177) 

(178) 
(179) 



(180) 



"Contrary to the spectral representation in previous sections, tlie limits of tlie integral are ±1/2 rather than ±00 since we are considering 
sampled (discrete) random fields. 
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On the other hand, using ^(x) = nis(x), we find that 



Moreover, it is easy to show that S{\) = —iF^{\) + iF (A), which yields 



rggW 



iV4 



E(l^^wi 



' + \FnW\' 



+ 0{N-^) 



(181) 



(182) 



The expectation of the error between the covariance of the infinite-length discrete-space and the periodic discrete 
Riesz transforms becomes 



E[Argg] ^ E[fggiO)] - E[fgg{0 



2 ^ /V2 ^ ^2 



N^ 



V-^l + ^^2 



2 \ sin^(7riV(A?ii - Ai)) sin^{TTN{Xn2 - Aa)) 
sin (7r(Ani — Ai)) sin (7r(An2 — A2)) 



5(A) dA. (183) 



Similar terms may be obtained for E[Ari^fi] and E[Argii]. 
Using the Taylor series of U given in ( |68| l, U2 is given by 

C/2- 



2n|Afgg -t- 2nf Af^,-^ - 4nin2Af~^ 



O (max(Ar3:^(:^) 



(184) 



Let us now write fff{0) = i? [?^//(0)] + dfff{0) and perform a Taylor series expansion of 1 + Sfff{0)/E [r//(0)] 
to find that 



U2 



2nlArgg + 2nf Af^^ - 4nin2A 



E[rffiO)] 
Now, taking the expectation of U2 it follows that 



E[rff{Q) 



O (max(Ar|'~^^~P 



^.,,1 2niii;[Af^g] + 2nf ii;[Af^^] - An,n2E[Afj^] ^ 



Plugging in the values for Afgg, . . ., it becomes 



E [U2] - 



1/2 



N/2-1 N/2-1 
^ ^ C(Ai,A2) 

Ai = -JV/2 A2=-A'/2 



S'(A)dA + C'(l/Af2), 



(185) 



(186) 



(187) 



E[fifmjo 

where C(Ai,A2) is defined in Lemma [T] with the exception that Ao is substituted by A. Now, using Lemma l2J it 
is possible to rewrite E [C/2] as 

E [U2\ = ^,^,. ,^,, / [sm^iTTC2)g{X, 711, 712) + sin2(^ci)g(A, 712, ni)] 5(A)dA + 0{1/N^). (188) 
NE[rff{0)]jQ 



Finally, taking the expectation also with respect to the direction, it becomes clear, using Lemmas pi ffl and Theorem 

[3] that the expectation is 

/■1/2 [41 4 1 41 

5'(A)dA 



E[U2 



4 14 1 4 
^A " 9 ^3 "^ 



N 



1/2 



0{1/N^ 



S{X)dX 



where we have used 



nl/2 



E[rff{0)]=rffiO)^2 f S{X)dX. 
Jo 

This is equivalent to ( |177| i because the method of moments estimator is unbiased for the variance. 



(189) 



(190) 
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